########利用optimize.leastsq函数进行拟合
import numpy as np 
from scipy import optimize
import matplotlib.pyplot as plt     #引入Excel库的xlrd
import sys

import xlrd                         #引入Excel库的xlrd

def GetXY(table, col, st_row, end_row=-1):
  "获取 x/Y 轴数组,从 col列 st_row行开始,默认到结尾"
  "数据为空则用前后数据中值替代"
  # return table.col_values(col-1)[st_row-1:end_row] #有空行没过滤
  lst = table.col_values(col-1)[st_row-1:end_row]
  for i in range(len(lst)):
    if(lst[i]==''):
      lst[i] = (lst[i-1]+lst[i+1])/2
  return lst

def leastsq(x,y):
  "最小二乘法拟合直线"
  "计算以p为参数的直线和原始数据之间的残差"
  def residuals(p): 
      k,b = p
      return y - (k*x + b)
  # l e a s t s q使得residuals()的输出数姐的平方和最小，参数的初始值为[1,0] 
  r = optimize.leastsq(residuals, [1, 0]) 
  k, b = r[0]
  print ("leastsq: k =" , k, "b =" , b)
  #####将计算结果可视化,画出拟合曲线和原数据点
  XX = np.arange(2250, 2900, 20)
  YY=k*XX+b
  plt.plot(XX,YY)
  plt.plot(x,y,'o')
  plt.xlabel('X')
  plt.ylabel('Y')
  # plt.show()



def line_fit_fun(x, k, b):  #直线拟合
    return k*x+b

def exp_fit_fun(x, a, b, c, dx): #指数拟合
    return a*np.exp(-b*(x-dx))+c

def log_fit_fun(x, a, b, c, dx): #log拟合
    return a*np.log(b*(x+dx))+c
  
def poly2_fit_fun(x, a, b, c, dx):  #多项式拟合
    return a*(x-dx)**2 + b*(x-dx) +c

def poly3_fit_fun(x, a, b, c, d, dx):  #多项式拟合
    return a*(x-dx)**3 + b*(x-dx)**2 +c*(x-dx) +d

def poly4_fit_fun(x, a, b, c, d, e, dx):  #多项式拟合
    return a*(x-dx)**4 + b*(x-dx)**3 + c*(x-dx)**2 + d*(x-dx) +e


def curve_fit(x, y, fit_fun, XX):
  "曲线拟合"
  "定义拟合函数"

  # #popt数组中，两个值分别是待求参数k,b
  popt, pcov = optimize.curve_fit(fit_fun, x, y, maxfev=10000000)

  # 根据不同拟合函数进行处理
  if( fit_fun == line_fit_fun ):
    print ("curve_fit：k =" , popt[0], "b =" , popt[1])
    YY = popt[0]*XX+popt[1]
    

  if( fit_fun == exp_fit_fun ):
    print ("curve_fit：a =" ,popt[0], "b =" ,popt[1], "c =" ,popt[2], "d =" ,popt[3])
    YY = popt[0]*np.exp(-popt[1]*(XX+popt[3]))+popt[2]
  if( fit_fun == log_fit_fun ):
    print ("curve_fit：a =" ,popt[0], "b =" ,popt[1], "c =" ,popt[2], "d =" ,popt[3])
    YY=popt[0]*np.log(popt[1]*(XX+popt[3]))+popt[2]
    
  if( fit_fun == poly2_fit_fun ):
    print ("curve_fit：k =" ,popt[0], "b =" ,popt[1])
    YY=popt[0]*(XX+popt[3])**2 + popt[1]*(XX+popt[3]) +popt[2]
    
  if( fit_fun == poly3_fit_fun ):
    for i in popt:
      print(i)
    print ("curve_fit：k =" ,popt[0], "b =" ,popt[1])
    YY=popt[0]*(XX+popt[4])**3 + popt[1]*(XX+popt[4])**2 +popt[2]*(XX+popt[4]) +popt[3]

  if( fit_fun == poly4_fit_fun ):
    for i in popt:
      print(i)
    print ("curve_fit：k =" ,popt[0], "b =" ,popt[1])
    YY=popt[0]*(XX+popt[5])**4 + popt[1]*(XX+popt[5])**3 +popt[2]*(XX+popt[5])**2 +popt[3]*(XX+popt[5]) +popt[4]

  # 计算方差
  mean = np.mean(y)  # 1.y mean
  ss_tot = np.sum((y - mean) ** 2)  # 2.total sum of squares
  ss_res = np.sum((y - fit_fun(x, *popt)) ** 2)  # 3.residual sum of squares
  r_squared = 1 - (ss_res / ss_tot)  # 4.r squared
  print("r_sq:", r_squared)


  #绘图
  plt.plot(x,y,'o')  
  plt.plot(XX,YY)

  plt.xlabel('X')
  plt.ylabel('Y')

  # plt.show()

def adc2flow( adc):
  p= [638.9,-1.074 ,0.000658,-1.763e-07, 1.772e-11]
  f=[]
  for x in adc:
    x2 = x*x
    x3 = x2*x
    x4 = x3*x
    f.append(p[4]*x4+p[3]*x3+p[2]*x2+p[1]*x+p[0])					
	
##    if(f<0): f = 0
##    elif(f>44): f = 44
  return f




if __name__ == '__main__':

##  adc = [i for i in range(2000,4000)]
##  flow = adc2flow(adc)
##  plt.plot(adc,flow)
##  plt.show()
  
  #导入需要读取Excel表格的路径
  # data = xlrd.open_workbook(r'.\flowadc-pwm.xlsx') #返回Book对象
  data = xlrd.open_workbook(r'.\flow-flowadc.xlsx') #返回Book对象
  table = data.sheets()[0]                      #读取第一张表

  flowAdc,flow = [],[]
  flowAdc = np.array(GetXY(table, col = 2, st_row = 2))#,end_row=62
  flow = np.array(GetXY(table,col=3,st_row=2))#,end_row=62

  # 拟合 流量-flowadc曲线
  print("flowadc-adc")
  curve_fit(flowAdc, flow, poly4_fit_fun, np.arange(2200, 3100, 0.1))
  plt.show()

  # 根据拟合曲线找到 流量0-25 对应大概的adc，ipython定义以下函数并找到对应adc
  def fun(x):
    a = -4.902515607926308e-11
    b = 5.411210946425652e-07
    c = -0.002211704816211876
    d = 3.9989549243865126
    e = -2708.0061578589834
    dx = 0.4863170267825967
    return a*(x-dx)**4 + b*(x-dx)**3 + c*(x-dx)**2 + d*(x-dx) + e


  # # 拟合 flowadc-流量曲线(不知为啥拟合有问题)
  # print("adc-flowadc")
  # curve_fit(flow,flowAdc,poly4_fit_fun,np.arange(0,60,0.1))



  # x=np.arange(2000,3400,2)
  # a=-101
  # b=-14
  # c=1400
  # dx=-2870
  # y=log_fit_fun(x,a,b,c,dx) #log拟合

  # for i in range(41):
  #   print(log_fit_fun(i,a,b,c,dx) )
  # # print(y)
  # plt.plot(x,y,'o')  
  
  plt.show()
